Monte Carlo simulace zdrojové přiměřenosti¶

připraveno pro výběrové řízení na pozici Analytik strategie ve společnosti ČEPS¶

Petr Louša¶

In [1]:
%matplotlib inline
In [2]:
import warnings

warnings.simplefilter(action="ignore", category=FutureWarning)

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go

from utils import *
In [3]:
pd.set_option("display.max_rows", 11)
pd.set_option("display.float_format", "{:.2f}".format)

Priprava dat¶

  • nactu data z Excelu
  • pripravim tabulku jednotlivych bloku
  • pridam casovou dimenzi
In [4]:
yearly_params, sources_definition = read_parameters("./parametry.xlsx", 3)
yearly_params = yearly_params.set_index("Rocni cyklus")
In [5]:
blocks = prepare_table_blocks(sources_definition)
blocks
Out[5]:
type block_no power priority shutdown_prob
0 FVE 0 2100.00 1 40
1 Jadro 0 716.67 2 10
2 Jadro 1 716.67 2 10
3 Jadro 2 716.67 2 10
4 Jadro 3 716.67 2 10
... ... ... ... ... ...
22 Plyn 3 287.50 4 20
23 Plyn 4 287.50 4 20
24 Plyn 5 287.50 4 20
25 Plyn 6 287.50 4 20
26 Plyn 7 287.50 4 20

27 rows × 5 columns

In [6]:
days = pd.DataFrame(pd.date_range(start="2022-01-01", end="2022-12-31", freq="D"), columns=["date"])
df_system = days.merge(blocks, how="cross")
df_system
Out[6]:
date type block_no power priority shutdown_prob
0 2022-01-01 FVE 0 2100.00 1 40
1 2022-01-01 Jadro 0 716.67 2 10
2 2022-01-01 Jadro 1 716.67 2 10
3 2022-01-01 Jadro 2 716.67 2 10
4 2022-01-01 Jadro 3 716.67 2 10
... ... ... ... ... ... ...
9850 2022-12-31 Plyn 3 287.50 4 20
9851 2022-12-31 Plyn 4 287.50 4 20
9852 2022-12-31 Plyn 5 287.50 4 20
9853 2022-12-31 Plyn 6 287.50 4 20
9854 2022-12-31 Plyn 7 287.50 4 20

9855 rows × 6 columns

Kalkulace dostupneho vykonu¶

  • vliv nahodnych odstavek
  • vliv pocasi na FVE
In [7]:
df_system = apply_shutdowns(df_system)
df_system
Out[7]:
date type block_no power priority shutdown_prob available_power
0 2022-01-01 FVE 0 2100.00 1 40 2100.00
1 2022-01-01 Jadro 0 716.67 2 10 716.67
2 2022-01-01 Jadro 1 716.67 2 10 716.67
3 2022-01-01 Jadro 2 716.67 2 10 716.67
4 2022-01-01 Jadro 3 716.67 2 10 716.67
... ... ... ... ... ... ... ...
9850 2022-12-31 Plyn 3 287.50 4 20 287.50
9851 2022-12-31 Plyn 4 287.50 4 20 287.50
9852 2022-12-31 Plyn 5 287.50 4 20 287.50
9853 2022-12-31 Plyn 6 287.50 4 20 287.50
9854 2022-12-31 Plyn 7 287.50 4 20 287.50

9855 rows × 7 columns

In [8]:
yearly_params
Out[8]:
leden cervenec zdroj
Rocni cyklus
Spotreba [GWh] 8500 6500 https://www.eru.cz/rocni-zprava-o-provozu-elek...
Utilizace FVE [%] 10 50 NaN
In [9]:
fve_utilization = yearly_params.loc["Utilizace FVE [%]"]

mask = df_system["type"] == "FVE"
df_system.loc[mask, "available_power"] = adjust_fve_power(
    df_system.loc[mask, "available_power"],
    january_utilization=fve_utilization["leden"],
    july_utilization=fve_utilization["cervenec"],
)

Spotreba¶

  • pridam denni spotrebu
  • vypnu nadbytecne zdroje
  • spocitam dny s nedodavkou
In [10]:
demand_params_MWh = yearly_params.loc["Spotreba [GWh]"] * 1000
demand = prepare_demand(days, demand_params_MWh["leden"], demand_params_MWh["cervenec"])
df_system = df_system.merge(demand, on="date")
df_system
Out[10]:
date type block_no power priority shutdown_prob available_power demand
0 2022-01-01 FVE 0 2100.00 1 40 210.00 11805.56
1 2022-01-01 Jadro 0 716.67 2 10 716.67 11805.56
2 2022-01-01 Jadro 1 716.67 2 10 716.67 11805.56
3 2022-01-01 Jadro 2 716.67 2 10 716.67 11805.56
4 2022-01-01 Jadro 3 716.67 2 10 716.67 11805.56
... ... ... ... ... ... ... ... ...
9850 2022-12-31 Plyn 3 287.50 4 20 287.50 11805.35
9851 2022-12-31 Plyn 4 287.50 4 20 287.50 11805.35
9852 2022-12-31 Plyn 5 287.50 4 20 287.50 11805.35
9853 2022-12-31 Plyn 6 287.50 4 20 287.50 11805.35
9854 2022-12-31 Plyn 7 287.50 4 20 287.50 11805.35

9855 rows × 8 columns

In [11]:
fig = plot_power(df_system, "available_power")
fig.show()

Vypnuti nadbytecnych zdroju¶

In [12]:
df_system.sort_values(by=["date", "priority", "available_power", "block_no"], inplace=True)
df_system["total_power"] = df_system.groupby("date")["available_power"].cumsum()

df_system["turn_off"] = (df_system["total_power"] - df_system["power"]) > df_system["demand"]
df_system = df_system.assign(final_power=lambda x: (~x.turn_off) * x.available_power)
df_system
Out[12]:
date type block_no power priority shutdown_prob available_power demand total_power turn_off final_power
0 2022-01-01 FVE 0 2100.00 1 40 210.00 11805.56 210.00 False 210.00
1 2022-01-01 Jadro 0 716.67 2 10 716.67 11805.56 926.67 False 716.67
2 2022-01-01 Jadro 1 716.67 2 10 716.67 11805.56 1643.33 False 716.67
3 2022-01-01 Jadro 2 716.67 2 10 716.67 11805.56 2360.00 False 716.67
4 2022-01-01 Jadro 3 716.67 2 10 716.67 11805.56 3076.67 False 716.67
... ... ... ... ... ... ... ... ... ... ... ...
9850 2022-12-31 Plyn 3 287.50 4 20 287.50 11805.35 11418.40 False 287.50
9851 2022-12-31 Plyn 4 287.50 4 20 287.50 11805.35 11705.90 False 287.50
9852 2022-12-31 Plyn 5 287.50 4 20 287.50 11805.35 11993.40 False 287.50
9853 2022-12-31 Plyn 6 287.50 4 20 287.50 11805.35 12280.90 True 0.00
9854 2022-12-31 Plyn 7 287.50 4 20 287.50 11805.35 12568.40 True 0.00

9855 rows × 11 columns

In [13]:
fig = plot_power(df_system, "final_power")
fig.show()

Nedodavka¶

In [14]:
df_missed = df_system.groupby("date")[["total_power", "demand"]].max()
df_missed = df_missed.assign(missed=lambda x: np.maximum(x.demand - x.total_power, 0))
df_missed.reset_index(inplace=True)
df_missed[df_missed["missed"] > 0]
Out[14]:
date total_power demand missed
0 2022-01-01 11510.00 11805.56 295.56
3 2022-01-04 11012.50 11803.70 791.20
9 2022-01-10 11148.36 11788.92 640.56
10 2022-01-11 10787.04 11785.03 997.99
12 2022-01-13 10802.26 11776.03 973.77
... ... ... ... ...
356 2022-12-23 9733.33 11788.92 2055.59
357 2022-12-24 9512.50 11792.41 2279.91
361 2022-12-28 11429.17 11802.26 373.10
362 2022-12-29 11510.56 11803.70 293.14
363 2022-12-30 11431.08 11804.73 373.65

63 rows × 4 columns

In [15]:
df_missed.loc[df_missed["missed"] > 0, "missed"].agg({"sum": "sum", "count": "count"})
Out[15]:
sum     59807.45
count      63.00
Name: missed, dtype: float64
In [16]:
# ax = df_missed.plot.line(x="date", y="missed")
# ax.set_ylabel("Missing power [MW]")
fig = px.bar(df_missed, x="date", y="missed")
fig.update_layout(
    title_text="Unsufficient supply",  # title of plot
    xaxis_title_text="Date",  # xaxis label
    yaxis_title_text="Missing power [MW]",  # yaxis label
)
fig.show()

Utilizace zdroju¶

In [17]:
df_utilization = (
    df_system.groupby("type")[["power", "available_power", "final_power"]]
    .sum()
    .assign(
        utilization_avail=lambda x: x.available_power / x.power * 100,
        utilization_final=lambda x: x.final_power / x.power * 100,
    )
)
df_utilization
Out[17]:
power available_power final_power utilization_avail utilization_final
type
FVE 766500.00 136191.73 136191.73 17.77 17.77
Jadro 1569500.00 1419000.00 1419000.00 90.41 90.41
Uhli 3431000.00 2266183.33 2057816.67 66.05 59.98
Plyn 839500.00 672750.00 216200.00 80.14 25.75
In [18]:
fig = plot_utilization(df_utilization)
fig.show()
In [19]:
df_total_energy = df_system.groupby("type")[["final_power"]].sum()
df_total_energy["final_energy_GWh"] = df_total_energy["final_power"] * 24 / 1000
df_total_energy["mix_percent"] = df_total_energy["final_power"] / df_total_energy["final_power"].sum() * 100
df_total_energy.reset_index(inplace=True)
df_total_energy
Out[19]:
type final_power final_energy_GWh mix_percent
0 FVE 136191.73 3268.60 3.56
1 Jadro 1419000.00 34056.00 37.06
2 Uhli 2057816.67 49387.60 53.74
3 Plyn 216200.00 5188.80 5.65
In [20]:
fig = plot_total_energy(df_total_energy)

Monte Carlo simulace¶

In [84]:
df_results_missed, df_results_power = run_monte_carlo("./parametry.xlsx", 30)

  0%|          | 0/30 [00:00<?, ?it/s]100%|██████████| 30/30 [00:00<00:00, 30.39it/s]
In [85]:
fig = px.scatter(df_results_missed[df_results_missed["missed"] > 0], x="date", y="missed", color="replica")
fig.update_layout(
    title_text="Energy not served",  # title of plot
    yaxis_title_text="Missing power [MW]",  # xaxis label
    xaxis_title_text="Date",  # yaxis label
)
fig.show()
In [86]:
fig = px.histogram(
    df_results_missed[df_results_missed["missed"] > 0].groupby("replica")["missed"].agg(["sum"]), x="sum", nbins=10
)
fig.update_layout(
    title_text="Energy not served", xaxis_title_text="Missing power [MW]", bargap=0.1  # title of plot  # xaxis label
)
In [87]:
fig = px.histogram(
    df_results_missed[df_results_missed["missed"] > 0].groupby("replica")["missed"].agg(["count"]), x="count", nbins=10
)
fig.update_layout(
    title_text="Days with missing power",  # title of plot
    xaxis_title_text="Number of days",  # xaxis label
    bargap=0.1,
)
In [88]:
fig = px.bar(
    df_results_power[["type", "final_energy_GWh"]]
    .groupby("type")
    .agg(["mean", "std"])
    .droplevel(level=0, axis=1)
    .reset_index(),
    x="type",
    y="mean",
    error_y="std",
    color="type",
)
fig.update_layout(
    title_text="Energy mix",
    yaxis_title_text="Total energy delivered [GWh]",
    xaxis_title_text="Source type",
    showlegend=False,
)
fig.show()